Multigrid Solution of a Path Integral Formulation for the Hydrogen Atomf 
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1. Introduction 

Most methods for Path Integral Monte Carlo (PIMC) simulations are inefficient for 
paths with large number of points. This is due mostly to the critical slow down (CSD) for 
lower dimensions and to the volume factor [5] for many dimensions. Furthermore, current 
applications of PIMC methods for solutions of Schrodinger equations for atomic systems is 
restricted because of the singularity of the potential functions. These systems are usually 
solved by Green's function Monte Carlo (GFMC) algorithms [2-3]. 

In this work (1) A. path integral formulation appropriate for Coulomb potentials is 
presented is Sec. 2. This formulation is based on a quadratic polynomial approximation 
of the classical motion. The resulting action integral can be expressed as a standard 
linear expression of the action integral with an effective potential. This effective potential 
depends on r 2 , where r is the mesh size (time step). (2) A multigrid algorithm using a 
"unigrid" approach is used for solving the hydrogen atom (Sec. 3). This algorithm uses 
a linear interpolation of changes to coordinates [5] , rather then constant interpolation [1] . 
The use of linear interpolation eliminates the CSD with simple V(l, 1) multigrid cycles. 
The integrated decorrelation time Tint is below 10 for all measured observables. In Sec. 4 
the numerical results with the multigrid algorithm are presented. For comparison, the 
same problem is solved with a staging algorithm [6-11]. 
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2. Path integral formulation for the hydrogen atom 

In this section r^, i = 0, . . . , N — 1 denotes the 3D coordinate vector (x i7 Zi) of the 
electron, r is the meshsize (timestep) of the path and N is the number of path points. 
Since only the ground state is considered the path is closed and = tq. Vi denotes the 
potential function at r$, and M is the mass of the electron. 

The standard first order approximation of the action integral for a path V is 

N-l N-l 

S 1 (F) = rJ2V l + (M/2r) £ (r m - r,) 2 . (2.1) 

i=0 i=0 

The relative probability for the path is, 

P(T)=e- s ^ r K 
However, this approximation fails for the Coulomb potential 

V(r) = -1/r. 

The reason is that when the Metropolis [4] algorithm is used to generate new paths, the 
path points tend to concentrate near r = 0, resulting in meaningless measurements of 
observables. This happens because near r = 0, the transition probability from large r is 
very close to 1. Similarly, since the transition probability from r = to large r is very 
small. Hence a path point can not make the transition to larger r, once trapped near 
r = 0. 

In terms of the classical path, the standard approximation (2.1) assumes that during 
the time r, the particle a) Move along the straight line connecting r^ and ri+i, and b) The 
time that the particle spends along any portion of that line is proportional to its length. 
It is mostly the second assumption that breaks down near the singularity, since classically 
if ri is near 0, the particle spends most of its time away from it. 

To make a better approximation of the classical motion, assumption a) above is pre- 
served, but instead of b) the motion is assumed to be in a constant force field. For a 
motion from r, to r^ + i the classical acceleration a,, is given by 

a, = -u(V(r i+1 ) - V(ri))/M\r i+1 - n\, (2.2) 

where u is the unit vector in the direction from r^ to r$+i. The position as a function of 
time is, 

r(t) = r, + (r i+ i - r,)((t - U)/t) + a^t - U){t -t t - r). 

This classical motion conserves energy. Therefore the total energy for the motion 
r^ — > r$+i can be taken as the energy at time £j, 

Ei = - 1/r, 

where is 

Vi = (r i+ i - Ti)/T - &iT. (2.3) 
Thus, the action integral for a path T is 

N-l 

S 2 (F)=rJ2E t (2.4) 

i=0 
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2.1 Simplified path integral formulation 

From (2.3), the expression for v 2 contains three terms: (i) (rj+i — ri) 2 /r 2 . This is 
the same expression for the squared velocity that appears in the standard approximation 
(2.1). (it) r 2 a 2 . From (2.2), in the limit r — > 0, this term for the Coulomb potential is 
0(r 2 /r 4 ). Similarly, in the same limit the term (m) — 2ai(ri + i — r^) is 0(l/r 2 ). 

Near r = 0, the term (it) is dominant. For large r, the standard term (i) is dominant. 
It is therefore reasonable to approximate the kinetic energy by preserving only terms (i) 
and (ii), and omitting (Hi) altogether. The action integral is then, 

N-l N-l N-l 

S(T) = (M/2r) £ (r i+ i - r,) 2 " r ^ + t(A/M) ^ r 2 /r 4 , (2.5) 

i=0 i=0 i=0 

where A is a positive constant. The last two terms can be regarded as a sum over an 
effective potential 

Vi(r, t) = -1/n + Ar 2 /Mrf. (2.6) 

Note that the last term in (2.6) acts as a repulsive potential, and thus prevents the 
path points from falling into r = 0. For a fixed u = At 2 /Mr(u) , r(u) decreases as yfr. 
With V the action integral is 

N-l N-l 

S(T) = (M/2r) (^+i " r i? + * E ^ ^ 

i=0 i=0 

This considerably simplifies the Monte Carlo sampling of paths, as its form is the same as 
the standard approximation (2.1). 
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3. Multigrid algorithm 

To eliminate the correlation between paths, a multigrid algorithm is used for gener- 
ating paths. The algorithm in this study implements a "unigrid" approach, because of 
its simplicity. In this approach there is only a single grid and coarse level operations are 
emulated on this grid. Such an algorithm in not very effective in reducing the volume 
factor (see. [5]) for a large number of dimensions. However, it is effective in eliminating 
the CSD (see Sec. 4), which is the most important factor in efficiency degradation for 
lower number of dimensions such as the hydrogen atom. For higher dimensions, a fully 
multigridded algorithm should be used. 

The finest level / = contains the whole set of points i = 0, . . . , N — 1 of the path, 
corresponding to times U = it. N is assumed to be an integral power of 2. A coarser level 
/ = 1, . . . , l c contains the set of ni = N/2 1 points i\ = n2 l , n = 0,1, ... ,ni, corresponding 
to the times t it = t it T. 

On the finest level, a relaxation sweep consists of changing the coordinates of each 
path point in turn, and accepting or rejecting the new position by the standard Metropolis 
algorithm [4]. 

To explain the relaxations on coarse grids, the following notation is introduced. Let 
^n,i2 be the set of path points corresponding to times tn,tn+i, . . . , U2. Let Sa,i2 be the 
action integral of Vn^: 

i2-l i2 

S a ,i2 = (M/2r) ]T (r i+1 - r,) 2 + r J] % (3.1) 

i=il i=il 

where V is given by Eq. (2.6). 

A relaxation sweep on level / > 0, is emulated by visiting all path points corresponding 
to that level. When the point i\ is visited, simultaneous changes of the coordinates are 
made to T il _ 2 i ^ l+ 2 l ■ These changes are given by 

S l = i = i l -2 l ,...,i l ,...,i l + 2 l (3.2) 

where Si is a constant, Si is the 3D vector (S x , S y , S z ) of changes to the space coordinates 
(x,y,z) and £ is the 3D vector (£ x , £ 2 ). Each of the components of £ is a random variable 
with uniform distribution in [—1, 1] (See also Sec. 4.1). If S a is the action of r i; _ 2 i ^+2' 
before changes to coordinates are made, and S b is the action after the changes, the new 
configuration is accepted according to the transition probability 

P ab = mm(l,e- (sb - sa) ). 

Paths are generated by cycling on all levels, usually with V(l,l) cycles. W-cycles or other 
cycles which increase the number of passes on coarse levels relative to fine levels (i.e. 
higher cycle indices, [5] works as well, and are effective in reducing the decorrelation times 
between paths. However, in the "unigrid" approach, the work spent on a coarse level is not 
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much less than the work on a fine level. Therefore there is not much to gain in efficiency 
by performing more relaxations on coarse levels. Note, however, that using higher cycle 
index with the "unigrid" approach is an effective research tool for studying the behavior 
of the algorithm. 

Measurements of observables can be done either on the finest or on coarser levels. 
Because of the high correlations between adjacent path points on the finest level, some 
work may be saved making measurements on coarser levels. However, observables that 
depend on correlations between adjacent path points, such as the the mean of squared 
velocity can not be measured on coarse levels. 

4. Numerical experiments 

Numerical results for the hydrogen atom with the multigrid algorithm are presented 
in Sec. 4.1. For comparison, results with a simplified staging algorithm ([9-11]) for the 
same problem are presented in (Sec. 4.2). 

All computations in this section are done with the effective potential (2.6). The results 
of the calculations are not very sensitive to value of A in (2.6). In all numerical experiments 
in this section A = 0.005 is always used. 
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4.1 Multigrid results 

In Table 1, the results with V(l,l) cycles with the "unigrid" approach (Sec. 3) are 
described for different number of levels. The meshsize of a level I with iVj points is r; = 
32/A z . 

In a relaxation sweep on the finest level (/ = 0) the (xi, yi,Zi) coordinates are changed 
simultaneously at each path point by adding S = (S x ,S y ,S z ) to the current values. S is 
calculated by S = £So- So is a constant scalar optimized to yield Metropolis acceptance 
ratio of ~ 0.5 on that level. On coarser levels, simultaneous changes are made to the 3 
space coordinates of more than one point, as described in Sec. 3. As with So, the scalar 
constant Si for / > is chosen so that the Metropolis acceptance ratio on that level is 
~ 0.5. 

The observables measured are the mean radius of the electron < r > , the mean inverse 
of the radius < 1/r > and the mean kinetic energy < e& >. < eu > is approximated by 
averaging the operator 



(M/2t 2 ) {{x i+1 - Xi)(xi - Xi_i) + (y i+1 - yi){yi - 2/*_i) + (z i+1 - Zi){zi - z*_i)) 



where r is the meshsize ([12]). For each observable the integrated decorrelation time {Tint) 
is listed. 



L N T 


< r > Tint 


< 1/r > Tint 


< efc > Tint 


6 128 0.2500 

7 256 0.12500 

8 512 0.062500 

9 1024 0.031250 
10 2048 0.015625 


1.431(7) 10.0 
1.459(6) 7.8 
1.489(7) 9.7 
1.493(4) 7.3 
1.498(3) 8.0 


1.130(4) 7.0 
1.064(4) 6.0 
1.019(3) 7.0 
1.009(2) 5.4 
1.004(3) 5.5 


0.80(1) 3.0 
0.89(1) 2.7 
0.92(1) 2.6 
0.99(1) 2.4 
0.98(1) 2.3 



L - number of levels 

A^o - number of points on finest level 

To - meshsize of finest level 

Table 1 - Multigrid V(l,l) cycles 



It is quite evident from the table that for each measured observable no significant CSD 
exists in the range of 128-2048 path points. 

In Table 2 a typical behavior of Si as a function of / is described. The results in 
the table are for the last case listed in Table 1 (L = 10, No = 2048). The Metropolis 
acceptance ratio (MAR) is listed in the last column. For fine levels, Si is proportional to 
the square root of the meshsize. For coarser levels with mean distance between grid points 
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comparable to < r >, 5i becomes smaller with increased meshsize. 



7 
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C 

Ol 


MAR 





0.125 


0.505 


1 


0.175 


0.509 


O 

A 




U.Olo 


3 


0.343 


0.516 


4 


0.480 


0.518 


5 


0.672 


0.511 


6 


0.941 


0.479 


7 


1.054 


0.487 


8 


0.922 


0.490 


9 


0.646 


0.515 



Table 2 - behavior of 5i 
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4.2 Staging algorithm results 

As a comparison with the performance of the multigrid algorithm the results of a 
staging algorithm for the same problem are described. This staging algorithm is similar to 
the algorithm used in [6] , [7] , and is a simplified version of the algorithm described in [8] . 

In a relaxation sweep, each point % of the path becomes an end point of a chain of 
adjacent p path points i, i + 1, . . . , % +p + l. The two end points % and % +p + 1 are held fixed 
and the coordinates of all the other chain points are changed simultaneously with means 
and variances of a Levy walk. The new chain configuration is then accepted or rejected. 
Few trials per chain may be repeated to increase efficiency. 

Table 3 describes the results of the staging algorithm for the hydrogen atom, with 
different number of path points (N), chain length (N c , including the end points) and trials 
per chain (N try )- The meshsize r is 32 /N. The last three columns lists 1) the compu- 
tational work per relaxation sweep (W r ). W r is proportional to the number of random 
numbers with uniform distribution in [0, 1] produced in a single relaxation sweep. 2) The 
Metropolis acceptance ratio (MAR) of chain states, and 3) The integrated decorrelation 
times (Ti n t) f° r measurements of < r >. 

As seen in Table 3, for large N, MAR is approximately 0.5 for N c = 5 and Nt ry = 2. 
Increasing these two numbers results only in very small decrease of Ti n t. Clearly, although 
the staging algorithm is much faster than a primitive Metropolis algorithm, it does not 
eliminate the CSD. 



N 


T 






W r 


MAR 


Tint <T> 


32 


1.0 


4 


1 


0.21 


0.32 


15.6 








2 


0.44 


0.67 


7.7 






5 


2 


0.64 


0.32 


7.8 








3 


0.95 


0.48 


7.7 






6 


4 


1.7 


0.24 


9.5 








5 


2.1 


0.32 


8.5 


64 


0.5 


4 


1 


0.44 


0.45 


33.0 






5 


2 


1.3 


0.46 


14.6 








3 


1.9 


0.67 








6 


3 


2.5 


0.31 










4 


3.3 


0.41 










5 


4.1 


0.50 


12.5 


128 


0.25 


4 


1 


0.86 


0.50 


120.0 






5 


1 


1.30 


0.25 










2 


2.5 


0.51 


75.0 


256 


0.125 


4 


1 


1.7 


0.51 








5 


1 


2.5 


0.27 










2 


5.1 


0.53 


350.0 


512 


0.0625 


5 


2 


10.0 


0.53 


1240.0 



Table 3 - Staging algorithm 
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